Визуализация и решение обыкновенных дифференциальных уравнений в R

library("knitr")
opts_chunk$set(
  # cache=FALSE,
               message=FALSE,warning=FALSE) 

library("ggplot2") # для построения графиков
library("rasterVis")
library("fields")
library("deSolve")
library("bvpSolve")

Пакет rasterVis предназначен для изображения данных на реальных географических картах, поэтому там нужно понятие проекции. Мы пока просто введем это шаманское заклинание

proj <- CRS('+proj=longlat +datum=WGS84')

Построим график векторного поля для системы:

\[ \left\{ \begin{array}{l} \dot{y}_1=y_2 \\ \dot{y}_2=y_1+\cos(y_2) \end{array} \right. \]

Задаем решетку, и рассчитываем \(\dot{y}_1\) и \(\dot{y}_2\) в точках решетки:

y1 <- seq(-6, 6, .05)
y2 <- seq(-6, 6, .05)
df <- expand.grid(y1=y1, y2=y2)
df$y1dot <- df$y2
df$y2dot <- df$y1+cos(df$y2)

Рассчитываем длины и углы для стрелочек, помещаем результат в объект Raster.

df$len <- sqrt(df$y1dot^2+df$y2dot^2)
df$angle <- atan2(df$y1dot,df$y2dot)

df2 <- df[c("y1","y2","len","angle")]

rast <- rasterFromXYZ(df2,crs=proj)

Строим классический график со стрелочками

vectorplot(rast,isField=TRUE)

Строим няку с капельками

streamplot(rast, isField=TRUE)

Простой график можно руками построить без доп. пакетов. При этом нам нужно самостоятельно уменьшить количество стрелочек.

y1 <- seq(-6, 6, .5)
y2 <- seq(-6, 6, .5)
df <- expand.grid(y1=y1, y2=y2)
df$y1dot <- df$y2
df$y2dot <- df$y1+cos(df$y2)
plot(df$y1,df$y2,pch=".")
arrow.plot( df$y1,df$y2,df$y1dot,df$y2dot,
            arrow.ex=0.03,length=0.05) 

Решим ОДУ с начальным условиями

Решим систему ОДУ с начальными условиями

Описываем саму систему:

eq1 <- function(t,y,parampampam) {
  return(list(c(
    y[2],
    y[1]+cos(y[2])    
  )))
}

Начальные условия:

y.start <- c(y1=1,y2=4) 

Точки, в которых компьютер будет считать функцию:

t <- seq(0,10,by=0.01)

Решаем

sol <- ode(y=y.start,times=t,func=eq1)
sol <- data.frame(sol)
head(sol)

##   time       y1       y2
## 1 0.00 1.000000 4.000000
## 2 0.01 1.040018 4.003678
## 3 0.02 1.080076 4.007785
## 4 0.03 1.120176 4.012326
## 5 0.04 1.160324 4.017305
## 6 0.05 1.200524 4.022725

qplot(data=sol,time,y1)

Функция ode возвращает матрицу, а для рисования графиков удобнее табличка с данными, data.frame. Строчка sol <- data.frame(sol) переделывает матрицу в таблицу с данными.

Решим систему ОДУ с краевыми условиями

Описываем саму систему:

eq1 <- function(t,y,parampampam) {
  return(list(c(
    y[2],
    y[1]+cos(y[2])    
  )))
}

Граничные условия:

y.start <- c(y1=1,y2=NA) 
y.final <- c(y1=42,y2=NA)

Точки, в которых компьютер будет считать функцию:

t <- seq(0,10,by=0.01)

Решаем

sol <- bvptwp(yini=y.start,yend=y.final,
           x=t,func=eq1,
           nmax=2000)
sol <- data.frame(sol)
head(sol)

##      x        y1        y2
## 1 0.00 1.0000000 -1.553150
## 2 0.01 0.9845193 -1.543001
## 3 0.02 0.9691398 -1.532904
## 4 0.03 0.9538610 -1.522860
## 5 0.04 0.9386824 -1.512868
## 6 0.05 0.9236035 -1.502928

qplot(data=sol,x,y1)

Бесплатное приложение. Изображение функций двух переменных

Есть несколько способов представить себе функцию от двух переменных, \(z(x,y)\):

• 3D график
• Линии уровня
• Векторное поле градиентов функции

Создаем data.frame с декартовым произведением двух векторов

df <- expand.grid(x=seq(-2, 2, .01), y=seq(-2, 2, .01))

Изобразим функцию \(z(x,y)=(3\cdot x^2+y)\cdot e^{-x^2-y^2}\).

Cоздаем переменную z как функцию от x и y

df$z <- with((3*x^2 + y)*exp(-x^2-y^2),data=df)

r <- rasterFromXYZ(df, crs=proj)

Линии уровня функции z

contour(r)

Капельки текущие по градиенту

streamplot(r)

Направление градиентов, заодно вид сбоку для графика функции

vectorplot(r)